POP Steering parameterization

Eddy diffusivity implementation by Bates, Tulloch, Marshall, and Ferrari

POP Steering Parameterization

Parameterizing the Eddy Length Scale

$$K=u_{rms}∗{\Gamma * \color{red}{L_{eddy}} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)} $$

The three length scales under consideration:

  1. Rossby Deformation Radius $\quad\quad L_r = {c_r \over |f|}$
  2. Rossby Equatorial Radius $ \quad\quad= L_{req} = {\sqrt{c_r \over 2\beta}}$
  3. Rhine's Scale $ \quad\quad L_{Rh} = {\sqrt{u_{rms} \over \beta}} \sim {\sigma_{vi} \over \beta}$

The Rossby Deformation Radius $L_r$ depends on the the Baroclinic wave speed $c_r$ and the Coriolis force $f$. The CESM value of the first Baroclinic wave speed is derived as per eq 2.2 in Chelton 1998.

First Baroclinic Wave Speed $c_r$ Chelton vs CESM

In [1]:
%pylab inline
%run steering_setup ;
basedir='/scratch/jet'
###nc = Dataset(basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.007.pop.h.nday1.yr.0249.nc')
nc = Dataset(basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.007.pop.h.nday1.0249-01.nc')
omega=7.2921e-5
deg2rad=pi/180
f=2*omega*sin(lat*deg2rad)
Populating the interactive namespace from numpy and matplotlib

/loan/jet/install/anaconda/lib/python2.7/site-packages/mpl_toolkits/__init__.py:2: UserWarning: Module argparse was already imported from /loan/jet/install/anaconda/lib/python2.7/argparse.pyc, but /loan/jet/install/anaconda/lib/python2.7/site-packages is being added to sys.path
  __import__('pkg_resources').declare_namespace(__name__)

In [2]:
### get Variable from file c_rossby is in cm/s convert to m/s
c_rossby = nc.variables['C_ROSSBY'][0,:,:]


#Resample (aka re-project, re-grid) the NCEP data to target grid. First with nearest neighbour resampling...
c_rossby_nearest = pyresample.kd_tree.resample_nearest(orig_def, c_rossby, \
        targ_def, radius_of_influence=500000, fill_value=None)
c_rossby_mps = c_rossby_nearest/100.

fig1=MapZoneContour(lon_targcyc,lat_targcyc,c_rossby_mps,figsize=(14,9),
                    title="CESM First Baroclinic Wave Speed (m/s)",
                    fmt='%.1f' )
c_r_chelton=Image(filename=basedir+'/steering/chelton_sfig1.jpg')
display(c_r_chelton)
/loan/jet/install/anaconda/lib/python2.7/site-packages/matplotlib/text.py:52: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal
  if rotation in ('horizontal', None):
/loan/jet/install/anaconda/lib/python2.7/site-packages/matplotlib/text.py:54: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal
  elif rotation == 'vertical':

First Baroclinic Rossby Radius Chelton vs CESM

CESM Rossby deformation radius = $L_{eddy}$

$L_r = {c_r \over |f|}$

$L_{req} = {\sqrt{c_r \over 2\beta}}$

$\require{cancel} \cancel {L_{Rh} = {\sqrt{u_{rms} \over \beta}} \sim {\sigma \over \beta}}$

$L_{eddy} = min (L_r, L_{req}, \cancel{L_{Rh}}).$

In [3]:
fcor = nc.variables['FCOR'][0,:,:]
fcor_nearest = pyresample.kd_tree.resample_nearest(orig_def, fcor, \
        targ_def, radius_of_influence=500000, fill_value=None)
absf=np.abs(fcor_nearest)
l_rossby_calc=c_rossby_nearest/absf

btp = nc.variables['BTP'][0,:,:]
btp_nearest = pyresample.kd_tree.resample_nearest(orig_def, btp, \
        targ_def, radius_of_influence=500000, fill_value=None)

l_rossbyeq_calc=np.sqrt(c_rossby_nearest/(2*btp_nearest))
l_eddy_calc=np.minimum(l_rossby_calc,l_rossbyeq_calc)

### convert regridded radius from cm to km
cm2km=1.e-5
l_eddy_calc_km=l_eddy_calc*cm2km
clevs=[10,20,30,40,50,60,80,100,150,230]
fig2=MapZoneContour(lon_targcyc,lat_targcyc,l_eddy_calc_km,figsize=(14,9),
                    levels=clevs,
                    title="CESM First Baroclinic Rossby Radius (km) $L_{eddy} \sim min(l_r,l_{req})$",
                    fmt='%.1i' )
l_r_chelton=Image(filename=basedir+'/steering/chelton_fig2.jpg')
display(l_r_chelton)
/loan/jet/install/anaconda/lib/python2.7/site-packages/numpy/ma/core.py:937: RuntimeWarning: overflow encountered in multiply
  result = self.f(da, db, *args, **kwargs)

Parameterizing the Zonal Eddy Phase Speed $\color{red}{c}$

$K=u_{rms}∗{\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - \color{red}{c}|^2 /u_{rms}^2 (z=0)} $

$\cancel{c = - \beta * {L_r^2}}$ $L_r$ too high at equator

$c = - \beta * L_{eddy}^2$

In [4]:
#beta
clevs=arange(1.e-14,25.e-14,1e-14)
fig,ax,cbar =MapContourg(lon_targcyc,lat_targcyc,btp_nearest,
                    levels1=clevs,
#                    norm=norm,
                    addzonal=True,
                    figsize=(14,5),
                    title="Beta (1/cm*s)")
l_eddy_calc_sq=l_eddy_calc*l_eddy_calc
clevs=np.logspace(1,15,40)
norm=LogNorm()
fig,ax,cbar =MapContourg(lon_targcyc,lat_targcyc,l_eddy_calc_sq,
                    addzonal=True,
                    levels1=clevs,
                    norm=norm,
                    figsize=(14,5),
                    title="$L_{eddy}^2 (cm^2)$ log scale")
#cbar.ax.get_xaxis().get_major_formatter().labelOnlyBase = True
#cbar.ax.get_yaxis().set_major_formatter(plt.LogFormatter(10,  labelOnlyBase=True))
#ax.yaxis.set_major_locator(ticker.MultipleLocator)
#ax.yaxis.set_major_formatter(LogFormatter())
#cbar.formatter.set_scientific(False)
#cbar.formatter.set_powerlimits((0, 0))
#cbar.update_ticks()
cbar.set_ticks(np.logspace(1,15,15))

c=btp_nearest*l_eddy_calc*l_eddy_calc # units=1/(cm s) * cm^2 = cm/s


clevs=arange(-10,20,1)
fig,ax,cbar = MapContourg(lon_targcyc,lat_targcyc,c,
                    addzonal=True,
                    levels1=clevs,
#                    norm=norm,
                    setover='darkred',
                    figsize=(14,5),
                    title="eddy phase speed c = beta*leddy**2 (cm/s) linear scale(NOTE I'm plotting positive c here - changed to negative after this plot")
cbar.set_ticks(np.linspace(-10,20,7))
c=-1.*c
/loan/jet/install/anaconda/lib/python2.7/site-packages/matplotlib/contour.py:1515: UserWarning: Log scale: values of z <= 0 have been masked
  warnings.warn('Log scale: values of z <= 0 have been masked')

Hughes Phase Speed (cm/s) from Tulloch Marshall Smith '09

In [5]:
hughes_c=Image(filename=basedir+'/steering/Hughes_phase_speed.png')
display(hughes_c)

Parameterizing the Zonal Velocity $\color{red}{u_{mean}}$

$K=u_{rms}∗{\Gamma * L_{eddy} \over (1 + b1 * |\color{red}{u_{mean}} - c|^2 /u_{rms}^2 (z=0)} $

In [6]:
#################  PLOT (u-c)^2  convert m^2/s^2 #######################
umean=nc.variables['U_MEAN'][0,:,:]
umean10_nearest = pyresample.kd_tree.resample_nearest(orig_def, umean,
                targ_def, radius_of_influence=500000, fill_value=None)

umean10_mps=umean10_nearest*1.e-2

clevs=np.linspace(-60,60,40)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,umean10_nearest,
                    addzonal=True,
                    figsize=(14,5),
                    levels1=clevs,
#                    setover='darkred',
#                    setunder='darkblue',
                    title="U zonal velocity (cm/s)")
#fig.subplots_adjust(left=0.15, right=.9, bottom=0.1, top=0.9);
cbar.set_ticks(arange(-60,60,10))
pylab.xlim([-20.,20.])


# get 2D versions of the lat and lon variables add longitude start point here!
lon_ecco_orig, lat_ecco_orig = np.meshgrid(lon_ecco[:], lat_ecco[:])

orig_ecco_def = pyresample.geometry.GridDefinition(lons=lon_ecco_orig, lats=lat_ecco_orig)

ecco_u_nearest = pyresample.kd_tree.resample_nearest(orig_ecco_def, ecco_u, \
        targ_def, radius_of_influence=500000, fill_value=None)
ecco_u_nearest_masked=np.ma.masked_where(ecco_u_nearest<-10, ecco_u_nearest, copy=True)

#print ecco_u_cyc_masked[:,0]

clevs=np.linspace(-60,60,40)
#ecco_u_nearest = pyresample.kd_tree.resample_nearest(orig_def, ecco_u, \
#        targ_def, radius_of_influence=500000, fill_value=None)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,ecco_u_nearest_masked*100,
                    addzonal=True,
                    figsize=(14,5),
                    levels1=clevs,
#                    extend='both',
                    title="ECCO U zonal velocity (cm/s)")
cbar.set_ticks(arange(-60,60,10))
pylab.xlim([-20.,20.])
Out[6]:
(-20.0, 20.0)
In [8]:
clevs=-np.logspace(3,-3,20)
norm=matplotlib.colors.SymLogNorm(linthresh=1e-4, vmin=-1000, vmax=-.001)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,c,
                    addzonal=True,
                    levels1=clevs,
                    norm=norm,
                    figsize=(14,5),
                    title="eddy phase speed c = beta*leddy**2 (cm/s) log scale")

# Customize y tick lables
yticks=[-1000,-500,-200,-100,-50,-20,-10,-5,-2,-1,-.500,-.200,-.1,-.05,-.02,-.01,-.005,-.002,-.001]
cbar.set_ticks(yticks)
cbar.set_ticklabels(yticks)

#ax.ticklabel_format(axis='y', style='sci', scilimits=(-2,2))
#majorFormatter = ticker.FormatStrFormatter('%.3f')
#majorLocator   = ticker.LogLocator(base=10.0, subs=[1.0], numdecs=6, numticks=15)
#ticker.MultipleLocator(20)
#cbar.ax.yaxis.set_major_locator(majorLocator)
#cbar.ax.yaxis.set_major_formatter(majorFormatter)
#ax.yaxis.get_major_formatter().set_powerlimits((0, 1))
pylab.xlim([-200.,0.])



uminusc=umean10_nearest-c  #units cm/s
clevs=np.linspace(0,25,25)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,uminusc,
                    addzonal=True,
#                    norm=norm,
                    levels1=clevs,
                    figsize=(14,5),
                    title="U-c (cm/s) (max = %i/min = %i) (c calculated using Leddy)"%(uminusc.max(),uminusc.min()))
yticks=linspace(0,25,6)
zonalticks=np.linspace(0,20,3)
cbar.set_ticks(yticks)
cbar.set_ticklabels(yticks)
pylab.xlim([0,25])
pylab.xticks(zonalticks)
Out[8]:
([<matplotlib.axis.XTick at 0x2b2223a77990>,
  <matplotlib.axis.XTick at 0x2b2222935990>,
  <matplotlib.axis.XTick at 0x2b224c086cd0>],
 <a list of 3 Text xticklabel objects>)
In [9]:
bates_uminusc=Image(filename=basedir+'/steering/bates2013_uminusc.png')
display(bates_uminusc)
In [10]:
uminusc_mps_sq=uminusc*uminusc*1.e-4   #units cm/s*cm/s*m/100cm*m/100cm
clevs=np.logspace(-4,1,40)
norm=LogNorm()
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,uminusc_mps_sq,
                    addzonal=True,
                    levels1=clevs,
                    figsize=(14,5),
                    norm=norm,
                    title=" ${(U-c)}^2 {(m^2/s^2)}$ max = %g/min = %g (log scale) "%(uminusc_mps_sq.max(),uminusc_mps_sq.min()))
#customize colorbar ticks
yticks=[1e-4,3e-4,1e-3,3e-3,1e-2,3e-2,1e-1,3e-1,1,3,10]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.06])
pylab.xticks([0,.03,.06])

clevs=np.linspace(0,.06,25)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,uminusc_mps_sq,
                    addzonal=True,
                    levels1=clevs,
                    figsize=(14,5),
                    title=" ${(U-c)}^2 {(m^2/s^2)}$ max = %g/min = %g "%(uminusc_mps_sq.max(),uminusc_mps_sq.min()))
#customize colorbar ticks
yticks=[0,.02,.04,.06]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.06])
pylab.xticks([0,.03,.06])
Out[10]:
([<matplotlib.axis.XTick at 0x2b224c564cd0>,
  <matplotlib.axis.XTick at 0x2b224c564990>,
  <matplotlib.axis.XTick at 0x2b2223f94210>],
 <a list of 3 Text xticklabel objects>)
In [11]:
bates_uminuscsq=Image(filename=basedir+'/steering/Bates_u_minus_c_sq.jpg')
display(bates_uminuscsq)

Parameterizing the root-mean-square eddy velocity $\color{red}{u_{rms}}$ and the Eady Growth Rate $\color{red}{\sigma_{vi}}$

$K=\color{red}{u_{rms}}∗{\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)} $

$u_{rms} = alpha*\color{red}{\sigma_{vi}}*L_{eddy}$

The original derivation of $\sigma_{vi}$:

$\sigma_{vi} = {f \over \sqrt{R_i}}$

with $R_i$ the vertically integrated (over the 100 – 2000 m depth range) Richardson number.

An alternate derivation of $\sigma_{vi}$ is:

$ R_i = {N^2 \over { ( \frac {\partial u}{\partial z} )^2 + ( \frac {\partial v}{\partial z} )^2} }$

$ N^2 = {{-g \over \rho_0 }\frac {\partial \rho}{\partial z}} $

After hydrostatic and geostrophic approximations

$f \frac {\partial v}{\partial z} = {{-g \over \rho_0 }\frac {\partial \rho}{\partial x}}; \quad\quad f \frac {\partial u}{\partial z} = {{g \over \rho_0 }\frac {\partial \rho}{\partial y}} $

so

$\frac {\partial v}{\partial z} = {{-1\over f}{g \over \rho_0 }\frac {\partial \rho}{\partial x}}; \quad\quad \frac {\partial u}{\partial z} = {{1\over f}{g \over \rho_0 }\frac {\partial \rho}{\partial y}} $

$\therefore$

$ R_i = {f^2 N^2 \over \underbrace{{ {g^2 \over \rho_0^2 } ( \frac{\partial \rho}{\partial y})^2 + {g^2 \over \rho_0^2 } ( \frac{\partial \rho}{\partial z})^2} }_\text{m^4}}\quad\quad = \quad\quad {f^2N^2 \over m^4} $

$\sigma_{vi} = {f \over \sqrt{R_i}}\quad = \quad {f \over \sqrt{{f^2N^2 \over m^4}}}\quad = \quad {{\cancel f m^2} \over \cancel f N}$

$RX_1 = RX_{east} = \Delta\rho_x = \rho_{i+1,j} - \rho_{i,j}$

$RY_1 = RY_{north} = \Delta\rho_y = \rho_{i,j+1} - \rho_{i,j}$

$RZ_1 = RZ_{k+1} = \Delta\rho_z = \rho_{k} - \rho_{k+1}$

$\displaystyle{1 \over L_{R_i}} \displaystyle\int_{2000m}^{100m} \left\lbrace { {-g\over\rho_0}{\frac {\partial \rho} {\partial z}} \over { {g^2 \over \rho_0^2 } \left[( \frac{\partial \rho}{\partial y})^2 + ( \frac{\partial \rho}{\partial z})^2\right] } } \right\rbrace dz$

Note: missing $f^2$ which will be cancelled when forming $\sigma_{vi}$

$\quad\quad$ so $\cdots$ this is not $R_i$

Implementation notes

For Level K

Numerator : Top $= -grav * RZ_{SAVE}(\cdots k+1) * dzwr(k)$

Denominator :

$ \begin{align} work1 = p25 & * ( RX(..,i_{east},k)^2 \\ & + RX(..,i_{west},k)^2 \\ & + RX(..,i_{east},k+1)^2 \\ & + RX(..,i_{west},k+1)^2 ) / DXT(i,j)^2 \\ \end{align} $

$ \begin{align} work2 = p25 & * ( RY(..,j_{north},k)^2 \\ & + RY(..,j_{south},k)^2 \\ & + RY(..,j_{north},k+1)^2 \\ & + RY(..,j_{south},k+1)^2 \\ \end{align} $

$work3 = {\left( TOP \over (grav^2*(work1+work2))\right)}*dzw(k)$

Notes:

1)Need to be careful at top and bottom of ocean
2)Accurate dzw(k) for each (i,j) to form $L_{R_i}$
3) When constructing $sigma$ itself, use $RZ_{SAVE}$ with a minimum N value
4) use eps2
In [12]:
sigma_old = nc.variables['SIGMA_AVG'][0,:,:]
sigma_old_nearest = pyresample.kd_tree.resample_nearest(orig_def, sigma_old, \
        targ_def, radius_of_influence=500000, fill_value=None)

sigma_avg1 = nc.variables['SIGMA_AVG1'][0,:,:]
sigma_avg1_nearest = pyresample.kd_tree.resample_nearest(orig_def, sigma_avg1, \
        targ_def, radius_of_influence=500000, fill_value=None)
clevs=arange(.005,.05,.0025)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,sigma_avg1_nearest*86400,
                levels1=clevs,
                figsize=(14,5),
                setunder='darkblue',
                setover='darkred',
                title="New Eady inverse time scale calc using alt sigma $(days^{-1}) \sim  {{m^2} \over N}$")
yticks=[.01,.02,.03,.04,.05]
cbar.set_ticks(yticks)

tulloch_eady=Image(filename=basedir+'/steering/tulloch_eady.png')
display(tulloch_eady)
In [13]:
#################  PLOT u_rms with Leddy #######################
cm2m=.01
alpha=5.
u_rms_leddy = alpha*sigma_avg1_nearest*l_eddy_calc    # units = 1*1/s*cm = cm/s
u_rms_leddy_mps=u_rms_leddy*cm2m

clevs=np.linspace(0,.4,40)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,u_rms_leddy_mps,
                    addzonal=True,
                    levels1=clevs,
                    figsize=(14,5),
                    title="$u_{rms} ({m \over s}) \quad alpha * \sigma * l_{eddy}$ (alpha=%i,max = %i/min = %i) "%(alpha,u_rms_leddy_mps.max(),u_rms_leddy_mps.min()))
#customize colorbar ticks
yticks=[0,.1,.2,.3,.4]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.2])
pylab.xticks([0,.1,.2])
Out[13]:
([<matplotlib.axis.XTick at 0x2b224cbc02d0>,
  <matplotlib.axis.XTick at 0x2b224cbc0a50>,
  <matplotlib.axis.XTick at 0x2b224ccc5750>],
 <a list of 3 Text xticklabel objects>)
In [14]:
bates_urms=Image(filename=basedir+'/steering/bates2013-urms.png')
display(bates_urms)
In [15]:
#################  PLOT u_rms^2  convert m^2/s^2 #######################
u_rms_leddy_mps_sq=u_rms_leddy_mps*u_rms_leddy_mps
clevs=np.linspace(0,.06,40)
norm=LogNorm()
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,u_rms_leddy_mps_sq,
                    addzonal=True,
                    levels1=clevs,
                    figsize=(14,5),
                    title="$u_{rms}^2 ({m^2 \over s^2})$ (alpha=%i,max = %g/min = %g) "%(alpha,u_rms_leddy_mps_sq.max(),u_rms_leddy_mps_sq.min()))
#customize colorbar ticks
yticks=[0,.02,.04,.06]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.06])
pylab.xticks([0,.03,.06])
Out[15]:
([<matplotlib.axis.XTick at 0x2b224c911650>,
  <matplotlib.axis.XTick at 0x2b224c911550>,
  <matplotlib.axis.XTick at 0x2b224c1a2cd0>],
 <a list of 3 Text xticklabel objects>)
In [16]:
urms_bates=Image(filename=basedir+'/steering/Bates_urms_sq.jpg')
display(urms_bates)
In [17]:
uminusc_sq_ovr_urms_sq = uminusc_mps_sq/u_rms_leddy_mps_sq             #units cm/s /  cm/s  
uminusc_sq_ovr_urms_sq_masked=np.ma.masked_where(uminusc_sq_ovr_urms_sq >100., uminusc_sq_ovr_urms_sq, copy=True)

clevs=np.logspace(-3,4,21)
norm=LogNorm()
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,uminusc_sq_ovr_urms_sq_masked,                
                    addzonal=True,
                    levels1=clevs,
                    norm=norm,
                    figsize=(14,5),
                    setover='darkred',
                    setunder='darkblue',
                    title="${(U-c)}^2 \over u_{rms}^2$ (masking values over 100)")                    
#customize colorbar ticks
yticks=np.logspace(-3,4,15)
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,100])
pylab.xticks([0,50,100])

clevs=arange(0,5,.1)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,uminusc_sq_ovr_urms_sq_masked,                
                    addzonal=True,
                    levels1=clevs,
                    figsize=(14,4),
                    setover='darkred',
                    setunder='darkblue',
                    title="${(U-c)}^2 \over u_{rms}^2$ (masking values over 100)")
#customize colorbar ticks
yticks=[0,2,4]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,8])
pylab.xticks([0,4,8])
Out[17]:
([<matplotlib.axis.XTick at 0x2b224e1f6f50>,
  <matplotlib.axis.XTick at 0x2b224e341710>,
  <matplotlib.axis.XTick at 0x2b224df50210>],
 <a list of 3 Text xticklabel objects>)
In [18]:
bates_uminusc_sq_ovr_urms_sq =Image(filename=basedir+'/steering/Bates_suppress.jpg')
display(bates_uminusc_sq_ovr_urms_sq )

Parameterizing the Mixing term $\color{red}{L_{mix}}$

$\color{red}{L_{mix}} = \Gamma * L_{eddy} * Suppression$

$\Gamma = 0.35$

$Suppression= {1 \over (1 + b1 * |u_{mean} - c|^2 /u_{rms (z=0)}^2 )}$

$\color{red}{L_{mix}} = {\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)}$

In [19]:
gamma=.35
b1=4.0
supp=1./(1+b1*uminusc_sq_ovr_urms_sq)
supp1b1=1./(1+1*uminusc_sq_ovr_urms_sq)

cm2m=.01
l_eddy_calc_m=l_eddy_calc*cm2m
lmix=gamma*l_eddy_calc_m*supp
lmix1b1=gamma*l_eddy_calc_m*supp1b1


clevs=np.logspace(0,5,40)
norm=LogNorm()
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,lmix,
                    addzonal=True,
                    norm=norm,
                    levels1=clevs,
                    figsize=(14,5),
                    title="$L_{mix} = {\Gamma * L_{eddy} * Suppression} $")

pylab.xlim([0,10e3])
cbar.set_ticks([1,3,1e1,3e1,1e2,3e2,1e3,3e3,1e4,3e4,1e5])



clevs=np.linspace(0,1,30)
norm=matplotlib.colors.Normalize()
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,supp,
                    addzonal=True,
                    levels1=clevs,
                    figsize=(14,5),
                    title="Suppression factor ${1 \over (1 + b1 * |u_{mean} - c|^2 /u_{rms (z=0)}^2)}$")

#customize colorbar ticks
yticks=np.linspace(0, 1, 6)
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.8])
pylab.xticks([0,.25,.5,.75])
Out[19]:
([<matplotlib.axis.XTick at 0x2b224f68e150>,
  <matplotlib.axis.XTick at 0x2b224f68e590>,
  <matplotlib.axis.XTick at 0x2b224f2a1e50>,
  <matplotlib.axis.XTick at 0x2b224f706a50>],
 <a list of 4 Text xticklabel objects>)
In [20]:
bates_suppression =Image(filename=basedir+'/steering/bates2013_suppression.png')
display(bates_suppression )

Parameterizing Eddy diffusivity $\color{red}{K}$

$\color{red}{K}=u_{rms}∗L_{mix}$

In [21]:
kdiff=u_rms_leddy_mps*lmix
kdiff1b1=u_rms_leddy_mps*lmix1b1

clevs=np.logspace(0,4,40)
norm=LogNorm()
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,kdiff,
                    addzonal=True,
                    norm=norm,
                    levels1=clevs,
                    figsize=(14,5),
                    title="Eddy diffusivity $K$ (log scale)")
#customize colorbar ticks
yticks=np.logspace(0, 4,5)
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,11e3])
pylab.xticks([0,5e3,10e3])

clevs=np.linspace(0,1.e4,40)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,kdiff,
                    addzonal=True,
                    levels1=clevs,
                    figsize=(14,5),
                    title="Eddy diffusivity $K$ ")
#customize colorbar ticks
#cbar.set_ticks([1,3,1e1,3e1,1e2,3e2,1e3,3e3,1e4,])
yticks=np.linspace(0, 10000,6)
print yticks
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,11e3])
pylab.xticks([0,5e3,10e3])


clevs=np.linspace(0,1.e4,40)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,kdiff1b1,
                    addzonal=True,
                    levels1=clevs,
                    figsize=(14,5),
                    title="Eddy diffusivity $K$ With a suppression scaling factor $b1$ of 1 instead of 4")
#customize colorbar ticks
#cbar.set_ticks([1,3,1e1,3e1,1e2,3e2,1e3,3e3,1e4,])
yticks=np.linspace(0, 10000,6)
print yticks
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,11e3])
pylab.xticks([0,5e3,10e3])
[     0.   2000.   4000.   6000.   8000.  10000.]
[     0.   2000.   4000.   6000.   8000.  10000.]

Out[21]:
([<matplotlib.axis.XTick at 0x2b224fd463d0>,
  <matplotlib.axis.XTick at 0x2b224fd46cd0>,
  <matplotlib.axis.XTick at 0x2b224f99b610>],
 <a list of 3 Text xticklabel objects>)
In [22]:
bates_k =Image(filename=basedir+'/steering/bates2013_K.png')
display(bates_k )

Steering Level status

1. L_eddy values ok
2. Rossby wave speed ok
3. Eady inverse timescale ok
    a. Use alternate calculation otherwise has problems at equator and
    b.Richardson averages are too high causing sigma to drop too low
4. Zonal Eddy Phase Speed (c) 
    a. Too high 10S:10N
    b. Missing westward velocities in ACC which appear 
        in Hughes Data used in this parameterization
5. (U-c) Because of problem with c this term is also too high 10S:10N.
6. U_rms  (alpha*sigma*l_eddy)  OK        
    a. Using alpha of 5
7. Zonal mean velocity looks ok compared to ECCO annual average
8. Suppression OK
9. K is low by a factor of 5

Todo

1. Need to work on Zonal phase speed c
2. Plot vertical distribution of K
In [23]:
urms_bates2=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-002.jpg')
urms_bates3=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-003.jpg')
urms_bates4=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-004.jpg')
urms_bates5=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-005.jpg')
urms_bates6=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-006.jpg')
urms_bates7=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-007.jpg')
urms_bates8=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-008.jpg')
urms_bates9=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-009.jpg')
urms_bates10=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-010.jpg')
urms_bates11=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-011.jpg')
display(urms_bates3,urms_bates4,urms_bates5,urms_bates6,urms_bates7,urms_bates8,urms_bates9,urms_bates10,urms_bates11)

CESM $U_{mean} = U_{resolved}$ integraged over {10,50,100,200,500} meters

In [24]:
umean_cesm=Image(filename=basedir+'/steering/umean.jpg')
display(umean_cesm)
In [25]:
bates_len_vel=Image(filename=basedir+'/steering/bates_len_vel.jpg')
display(bates_len_vel)